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SECTION  I 
INTRODUCTION 

The  accurate  solution  of  fluid  dynamics  problems  of  practical 
significance  continues  to  comprise  a challenging  area  of  research 
owing  to  the  complexity  of  the  differential  equations  and  boundary 
conditions  encountered  in  the  matliematical  representation  of 
these  problems.  At  the  same  time,  the  demands  placed  on  technology 
today  require  attacking  problems  of  increasing  complexity.  Rapid 
advances  in  computing  machines  and  the  progress  in  the  development 
of  numerical  algorithms  now  allow  the  consideration  of  some 
complex  problems  not  heretofore  solvable.  However,  the  formulation 
of  the  problem  itself  has  a very  significant  role  in  the  success 
of  its  solution  and  in  the  efficiency  with  which  the  solution  is 
obtainable.  The  emphasis  o\  efficiency  of  a solution  increases 
with  the  complexity  of  the  problem  as  it  can  frequently  make  the 
difference  between  the  practicality,  or  otherwise,  of  obtaining 
the  solution  of  the  problem. 

An  essential  step  in  the  proper  formulation  of  the  problem 
consists  of  the  selection  of  an  appropriate  system  of  coordinates 
for  representing  the  problem.  In  order  for  a coordinate  system 
to  be  suitable,  it  must  satisfy  several  fundamental  requirements. 
Primary  among  these  is  the  requirement  that  the  coordinates  be 
surface-oriented  along  all  the  boundaries  of  the  problem,  as  this 
leads  to  simple  and  accurate  expressions  for  the  boundary 
conditions  of  the  problem.  For  flow  configurations  with  boundaries 
having  analytical  shapes,  surface-oriented  coordinates  may  be 


determined  through  analytical  transformations.  The  technique  of 
conformal  transformations  has  been  recently  used  very  ingeniously 


to  formulate  surface-oriented  coordinates  for  some  non-simple 

(1)  * 

geometries.  The  procedure  developed  by  Ives  using  a sequence 

of  conformal  transformations  was  later  employed  by  Grossman  and 
■ (2) 

Melnik'  ' to  successfully  compute  the  transonic  flow  over  two- 

13) 

element  airfoil  systems.  Moretti  has  developed  a transformation 
technique,  consisting  of  repeated  applications  of  a simple 
conformal  mapping  functioxi,  to  generate  surface-oriented  coor- 
dinates for  complex  geometries,  such  as  a typical  airframe 
configuration. 

Of  course,  it  is  always  desirable  to  seek  analytical  trans- 
formation techniques  rather  than  numerical  transformations 
because,  in  the  case  of  the  former,  the  mapping  as  well  as  all 
coordinate-derivatives  are  defined  exactly.  Nevertheless, 
analytical  transformations  have  certain  inherent  limitations. 

The  conformal  mapping  technique,  which  has  been  so  successful  for 
two-dimensional  transformations,  is  not  defined  for  three- 
dimensional  transformations.  It  can,  however  be  used  for  a 
three-dimensional  problem  to  determine  surface-oriented  two- 
dimensional  coordinates  in  planes  perpendicular  to  the  third 
dimension  in  the  problem.  This  approach  is  useful  when  the 
variations  of  the  body  geometry  in  this  third  dimension  are  only 
gradual.  Thus,  it  is  clear  that,  for  general  three-dimensional 
configurations  as  well  as  for  arbitrary  two-dimensional  configura- 
tions, the  more  general  numorical  transformation  techniques  need 


* Superscript  numerals  denote  similarly  numbered  references  at 
end  of  report. 
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to  be  investigated. 

Considerable  progress  has  been  made  in  the  development  of 
numerical  transformations  for  determining  surface-oriented 
coordinates  for  complex  configurations  [e.g.,  Refs.  4 and  5,  as 
well  as  the  additional  references  listed  therein] . While  the 
fundamental  procedure  of  these  numerical  transformations  has 
been  well-known  for  more  than  a decade,  it  has  been  given  wide 
exposition  only  recently.  Thompson  et  al.  employed  the 
technique  to  formulate  the  coordinates  for  several  external-flow 
configurations,  involving  doubly  connected  regions.  Ghia 
et  al.  ^ developed  the  analysis  for  three-dimensional  numerical 
transformations,  with  the  cylindrical  coordinate  system  as  the 
starting  system;  for  geometrical  similarity  with  respect  to  the 
radial  direction,  the  coordinate- transformation  equations  can  be 
made  to  resemble  those  for  the  two-dimensional  case. 

The  procedures  developed  are  applicable  to  arbitrary 
geometries  and,  for  two-dimensional  configurations,  consist,  in 
principle,  of  solving  elliptic  equations  of  the  form 

V^n  = Q , = R , (1) 

after  expressing  them  in  their  inverted  forms,  i.e.,  with  the 
roles  of  the  independent  and  the  dependent  variables  interchanged. 
The  boundary  conditions  needed  for  a unique  solution  of  Equation 
(1)  consist  of  specified  distributions  of  mesh  points  along  all 
boundaries  of  the  solution  domain..  In  Equation  (1)  , the  symbol  n 
denotes  the  transverse  coordinate,  and  c denotes  the  coordinate 
along  the  surface  [Figure  1] . The  forcing  functions  Q and  R 
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FIG.  1,  SCHEMATIC  REPRESENTATION  OF  ARBITRARY  GEOMETRICAL 
CONFIGURATION  - PHYSICAL  PLANE 
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serve  to  control  the  coordinate  spacing  in  the  interior  of  the 
physical  domain.  The  resulting  solution  can  provide  a more 

1 

optimum  distribution  of  a given  number  of  mesh  points  than  would 
be  possible  v.’i'.h  a conformal  coordinate  system  using 
* Q = R = 0.  In  the  computational  domain,  the  transformed  coor- 

j dinates  are  rectangular  and  uniformly  spaced. 

While  the  basic  procedure  is  quite  well  understood,  the 
transformation  involves  several  areas  of  detail  that  need  to 
p be  refined  and  formalized.  For  example,  considerations  of 

accuracy  of  the  flow  solution  require  that  regions  of 

! 

large  flow  gradients  be  appropriat ^ly  stretched  in  the  compu- 
tational domain.  This  is  to  be  achieved  by  employing  the  proper 

i 

j ' form,  in  magnitude  and  sign,  for  the  forcing  functions  in  the 

I 

coordinate  equations.  Although  it  is  simple  to  establish  the 
signs  of  the  forcing  functions,  the  magnitudes  usel  for  the  forcing 
functions,  in  order  to  achieve  a desired  mesh-point  distribution, 
are  frequently  determined  by  a process  of  trial  and  error.  It  is 
necessary  to  develop  a systematic  procedure  for  determining  the 
forcing  function  needed  to  provide  the  required  spacing  of  the 
coordinate  lines.  Moreover,  when  a nonuniform  grid  is  desired 
{ in  the  physical  domain,  the  distribution  of  points  along  the 

boundaries,  prescribed  as  boundary  conditions,  must  be  consistent 
with  the  forcing  functions  used  for  the  solution  in  the  interior. 
Failure  to  maintain  this  consistency  may  lead  to  coordinates  that 
exhibit  discontinuities  in  their  derivatives  near  the  boundaries. 

Another  area  requiring  further  study  deals  with  the  question 
of  the  location  of  the  far-field  boundary  for  subsonic  external 
flow  past  a body.  One  answer  to  this  question  consists  of  locating 

k 

■ -A: 
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the  far-field  boundary  infinitely  far  from  the  body.  However, 
this  necessitates  the  use  of  an  analytical  transformation  to  map 
the  unbounded  physical  domain  to  a finite  region.  To  the  best 
knowledge  of  the  authors,  this  has  not  been  done  so  far  for 
arbitrary  bodies,  in  conjunction  with  a nximerical  transformation. 
Frequently,  the  approach  used  for  locating  the  far-field 
boundary  for  flow  past  a body  has  consisted  of  repeating  the 
calculations  with  two  or  more  different  far- field-boundary 
locations  until  the  resulting  differences  in  the  normal  or 
tangential  stress  distributions  on  the  body  surface  are  within  a 
certain  prescribed  limit.  As  discussed  later  in  the  present 
study,  this  process  may  not  be  necessary,  especially  for  the 
commonly  studied  bodies  such  as  circular  cylinders,  airfoils,  etc. 

It  is  also  necessary  to  recall  that  the  basic  approach  for 
numerical  generation  of  surface-oriented  coordinates  does  nothing 
to  maintain  or  ensure  orthogonality  for  the  new  coordinates. 
However,  the  use  of  an  orthogonal  or  near-orthogonal  coordinate 
system  is  often  advantageous  as  regards  convergence  rate  of 
the  numerical  scheme.  An  orthogonal  system  is  also  desirable 
if  one  intends  to  make  a direction-related  approximation  in 
the  flow  equations,  e.g.,  parabolization  with  respect  to  the 
strecimwise  direction,  at  least  in  the  viscous  region  near  the 
body  surface.  Therefore,  it  would  be  desirable  to  develop  a 
procedure  for  orthogonalization  of  a given  numerical  non-orthogonal 
coordinate  system.  This  may  be  included  as  an  optional  step  at 
the  end  of  the  main  solution  procedure. 

The  purpose  of  the  present  study  is  to  optimize  the  nuiiierical 
transformation  procedure  for  determining  a suitable  system  of 
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surface-oriented  coordinates  for  arbitrary  geometrical  con- 
figurations. The  resulting  coordinates  are  required  to  provide 
appropriate  stretching  of  the  viscous  region  near  the  body  surface, 
in  the  computational  domain,  so  as  to  be  suitable  for  high- 
Reynolds-niamber  flow  calculations.  It  is  also  desired  that  the 
coordinates  be  nearly  orthogonal  in  the  physical  domain,  at 
least  near  the  boundaries.  Attempt  is  made  to  provide  as  much 
formalism  to  the  procedure  as  possible,  so  as  to  maintain 
versatility  and  generality  for  application  and  minimize  the 
need  for  a trial-and-error  approach.  At  the  stime  time,  this 
should  not  be  interpreted  as  de-emphasizing  the  significance  of 
good  engineering  intuition  in  the  efficient  solution  of  a 
practical  problem.  In  addition  to  analyzing  those  aspects  of 
the  transformation  discussed  in  the  preceding  paragraphs,  effort 
is  also  made  to  improve  the  numerical  scheme  itself.  The  major 
steps  used  to  achieve  the  objectives  of  this  research  are 
discussed  in  the  following  sections  of  this  report. 


7 


SECTION  II 


DESCRIPTION  OF  BASIC  PROCEDURE  P^OR  NUMERICAL 
GENERATION  OF  SURFACE-ORIENTED  COORDINATES 


In  order  to  provide  the  proper  perspective  for  the  optimi- 
zation steps  to  be  discussed  in  this  report,  a brief  description 
is  first  given  of  the  basic  procedure  for  numerical  generation 
of  surface-oriented  coordinates  for  arbitrary  geometries.  The 
analysis  presented  here  is  for  two-dimensional  geometries  or  tor 
cylindrical  geometries  with  radial  similarity.  Therefore,  in 

the  physical  domain,  the  configuration  is  represented  in  terms 
of  the.  coordinates  (<ti,  z)  such  that 

{y  for  two-dimensional  configurations 

r0  for  cylindrical  configurations. 

Nevertheless,  all  steps  in  the  analysis  are  either  applicable 
or  extendable  to  three-dimensional  transformations. 

In  principle,  the  nijmerical  transformation  procedure  con- 
sists of  determining  the  boundary-oriented  transformed  coordinates 
as  the  solutions  of  the  following  Poisson  equations: 


1 

and 

C 

where 


+ "zz  - 0 <2> 

+ 'izz  - (3) 

n,  C are  the  transformed  surface-oriented  coordinates; 
iji , z are  the  physical  coordinates; 

Q,  R are  continuous  functions  of  (n,  C)  or  (iji,  z)  , 
specified  so  as  to  provide  desired  control  on 
the  coordinate  distribution. 
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The  boundary  conditions  for  Equations  (2)  and  (3)  must  be 
specified  on  the  basis  of  the  particular  geometry  under  con- 
sideration. For  non-simple  geometries,  they  will  generally  not 
be  along  constant  values  of  <|>  or  z. 

Therefore,  it  is  necessary  to  derive  a set  of  equations 
equivalent  to  Equations  (2)  and  (3),  such  that  only  n,  C 
appear  as  independent  variables.  This  has  been  done  previously 
and  results  in  the  following  'inverted'  equations: 


and 


a <1.^^  + 2b  + c + J'^(Q  + R = 0 


az  +2bz^+c  z^^  + J (Q  z + R z,)  “=  0 

nn  n?  CC  n C 


(4) 


(5) 


where 


2 2 
a = (t  + z 


b = + z^z^) 


C = (()2  + z^ 

n n 


(6) 


and  J=(|)z^-«})^z 
n C n 


(7) 


It  is  noted  that  the  quantity  J defined  by  Eq.  (7)  is  exactly 
the  Jacobian  of  the  coordinate-transformation. 

The  boundary  conditions  for  Equations  (4)  and  (5)  can  now 
be  specified  as  the  prescribed  values  of  (j)  and  z,  along  the 
boundaries  11=11,,  ii„  and  C=Ci  , Co  ; [see  Figure  2].  However, 
when  a non-uniform  spacing  of  the  transformed  coordinates  is 


desired  in  the  physical  plane,  the  boundary-values  of  41  and  z 


may  no  longer  be  prescribed  independent  of  the  forcing  functions 
Q and  R used  in  Eqns.  (4)  and  (5).  Until  now,  the  situation  was 
handled  by  giving  up  the  Dirichlat  boundary  conditions  for 
Eqns.  (4)- (5),  in  favor  of  Neumann  boundary  conditions  consisting 
of  zero  normal  gradients  of  4>  and  z at  the  boundaries.  The 
drawback  of  this  approach,  as  well  as  the  present  improved 
approach  for  treating  this  and  several  other  aspects  of  the 
general  transformation  procedure,  are  discussed  next.  It  is 
important  to  note  that  improvement  in  the  analysis  also  leads, 
frequently,  to  improvements  in  the  finite-difference  method 
for  the  numerical  solution  of  the  transformation  equations. 
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SECTION  III 


OPTIMIZATION  OF  NUMERICAL  COORDINATE- 
TRANSFORMATION  PROCEDURE 

1.  REPRESENTATION  OP  VISCOUS  REGION  NEAP  BODY  OR  OTHER  HIGH- 
G1^.DIENT  REGIONS 

As  the  Reynolds  numlaer  Re  for  a flow  problem  increases,  the 
thickness  of  the  viscous  region  near  the  body  surface  decreases, 
so  that  the  flow  variables  exhibit  large  transverse  gradiei^ts 
in  this  important  region.  Therefore  an  accurate  computation  of 
the  flow  solution  requires  that  this  region  be  appropriately 
stretched  in  the  transformed  domain  so  as  to  provide  a reasonable 
number  of  computational  points  within  this  region.  As  mentioned 
in  Section  1,  this  is  achieved  by  using  the  proper  form  for 
the  function  Q in  the  coordinate-transformation  equations.  The 
proper  sign  for  Q can  generally  be  determined  on  the  basis  of  the 
theory  of  harmonic  functions.  However,  no  formal  theory  is 
available  for  determining  the  magnitude  of  Q.  It  is  only  known 
that  the  magnitude  required  for  Q generally  increases  with  the 
degree  of  stretching  required  and,  hence,  with  increasing  Re.  A 
procedure  has  been  developed,  in  the  present  work,  for  determining 
Q(n)  such  that  a prescribed  nvunber  of  computational  points  are 
located  in  a Blasius  boundary  layer  at  a body  surface.  Before 
describing  this  procedure,  it  is  necessary  to  consider  some  of 
the  adverse  effects  of  the  use  of  a non-zero  Q on  the  computational 
efficiency  of  the  solution  of  the  courdinate-trans formation  equa- 
tions. As  the  following  analysis  and  discussion  shows,  most  of 
these  adverse  effects  have  been  minimized  in  the  present  work.  In 
fact,  the  analysis  has  shed  light  on  some  very  useful  aspects  of 
the  probler.. 
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a.  Boundary  Conditions  for  Coordinate  Transformation  with 


Forcing  Function 

It  had  been  observed  that  the  convergence  rate  of  the 
coordinate  solution  was  usually  retarded  by  the  use  of  a non- 
zero Q,  even  if  only  moderate  magnitudes  were  used  for  Q.  One 
of  the  causes  for  this  seems  to  b(j  the  following.  In  order  to 
avoid  discontinuities  in  the  transformed-coordinate  derivatives, 
it  is  necessary  that  the  distribution  of  points  along  the  end 
boundaries  AB  and  CD  [Fig.  1]  be  consistent  with  the  function 
Q used  in  the  interior.  The  method  used  so  far  to  achieve  this 
consisted  of  specifying  conditions  only  on  the  coordinate 
derivatives  rather  than  on  the  locations  of  the  points  along 
AB  and  CD,  This  leads  to  slower  convergence  rate  for  the 
solution,  especially  if  an  explicit  numerical  scheme  is  employed. 
Moreover,  when  derivative  boundary  conditions  were  used  at  both 
boundaries  ^B  and  CD,  the  solution  was  sometimes  reported®  to 
violate  the  maximum  principle,  resulting  in  cross-over  of  the 
coordinates. 

In  the  present  work,  use  is  made  of  the  limiting  form  of 
the  coordinate  equations at  the  end-boundary  (e.g.,  A3)  to 
determine,  prior  to  the  complete  solution,  the  point-distributioti 
at  this  boundary,  consistent  with  the  Q used  in  the  interior, 

Tne  end  boundary  is  generally  a straight  line,  though  not  neces- 
sarily a vertical  straight  line.  However,  the  required  point- 
distribution  may  be  determined  for  a vertical  line  and  then 
transformed  to  yield  the  distribution  along  a general  line 
such  as  A' B' . 
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FIG.  3.  TREATMENT  OF  BOUNDARY  CONDITIONS  ALONG  A STRAIGHT 
BOUNDARY. 

Along  the  boundary  AB,  where  t=constant,  the  coordinate  equa- 
tions (1)  reduce  to 

= Q(n)  (8) 

if  the  (n-C)  system  is  required  to  be  orthogonal  along  AB. 

Since  Q will  not,  in  general,  be  a simple  linear  function  of  rif 
Eqn.  (8)  cannot  be  integrated  directly.  Howesver,  in  its 
inverted  form,  this  equation  becomes 
3 


nn  n 


= 0 


(9) 


with  the  boundary  conditions 

4 (0)  = 0 ; 4 (1)  = 4 


max 


(1C) 


Equation  (9)  is  a nonlinear  differential  equation  for  4,  with  a 
varying  coefficient,  and  can  be  easily  solved  numerically,  using 
the  boundary  conditions  (10).  A simple  tridiagonal  scheme,  with 


linearization  or  quasi-linearization  of  the  term  4 , was  used 


to  solve  this  equation,  with  the  nonlinearity  being  restored 
through  an  iterative  procedure.  Figure  4 shows  the  solution, 
hence  the  point-distribution  along  AB,  obtained  with  Q(ri) 
given  by  the  following  analytical  form; 


® I ZT  (2bJ)  ] (11) 

k=l  ^ 

where  a^  < 0 , ] a^ 1 = a2  , 

= b^  = 0.1 

= 0 , 

TI2  = 1 . (12) 

This  form  for  Q(ri)  provides  a fine  mesh  near  as  well  as 
near  n2* 

The  resulting  distribution  is  then  transformed  to  a corre- 
sponding distribution  along  line  A'B'  such  that  incremental 
distances  along  A'B'  are  the  same  as  the  corresponding  increments 
along  AB.  Thus,  between  the  two  points  i and  (i+l) , 


((})..,  - (j)  . ) ^ + (Z'.T  ■■  Z*)^ 
'^1+1  ' 1+1  1 


(13) 


where  denotes  the  solution  of  Equation  (9)  along  AB,  Denoting 
the  slope  of  line  A'B'  by  m.  Equation  (13)  can  be  rewritten  as 

' (1  + ,14, 


where  m 


(4>  — )/  (z 

'’^i+l  ^i'^'  i+l 


z. ) . 
1 
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LUTION 


Equation  (14)  yields  the  z-coordinates  of  the  required  point- 
distribution  along  the  actual  end-boundary  A'B'  as 


z 


i+1 


(15) 


where  the  sign  to  be  used  depends  upon  the  sign  of  m.  If 
m > 0,  the  positive  sign  must  be  used,  while  the  negative  sign 
is  needed  if  m < 0.  Therefore,  Eq.  (15)  must  be  represented 
as  follows: 


z 


i+1 


z . 
1 


+ { 


(^i+1  - *i^ 

A+m^ 


(16) 


Thereafter,  may  be  determined  from  the  slope  m as 

({>  = <J) . + m(z  - z , ) . (17) 

1+1  1 1+1  1 


Thus,  starting  with  the  coordinates  of  A' , the  value  of 
({i  and  z for  the  required  points  along  A'B'  can  be  computed  using 
the  solution  of  Equation  (9)  and  the  relations  (16)- (17).  It 
must  be  noted  that,  if  m = 0,  Eqns.  (16)  and  (17)  must  be  replaced 

_ _ 

"i+i " "i  '*1+1  • *1*  ' 


♦■X,  ' 

1+1  1 


(18) 


Also,  if  m ->■  «,  the  proper  relations  are 


z . , 1 = z . 
1+1  1 


=41.  + (4’..,  “<!'.) 
1+1  1 1+1  1 


(19) 
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It  is  important  that  the  length  of  AB  along  which  Eqn.  (9) 
is  solved  must  be  exactly  equal  to  the  length  of  the  actual 
boundary  A'B'.  If  this  is  not  maintained,  the  nonlinearity 
of  Eqn.  (9)  leads  to  a point-distribution  along  A’B'  that 
corresponds  to  an  effective  function  Q' , whose  values  are  related 
to  Q by  the  square  of  ratio  of  the  lengths  of  AB  and  A'B'. 

b . Determination  of  Forcing  Function  from  Specified 
Dirichlet  Boundary  Conditions 

In  another  situation,  the  distribution  of  the  points  along 
the  boundaries  of  an  arbitrary  geometrical  configuration  may  be 
prescribed  on  the  basis  of  some  independent  physical  considera- 
tions for  the  flow  problem.  In  order  to  obtain  the  transformed 
coordinates  satisfying  these  boundary  conditions,  it  is  now 
necessary  to  determine  the  forcing  function,  which  is  con- 
sistent with  the  specified  boundary-point  distribution.  In 
doing  so,  it  is  also  desirable  to  maintain  local  orthogonality 
for  the  coordinates.  Though  applicable  in  general,  the 
procedure  described  below  discusses  the  determination  of  the 
forcing  function  R in  Eqn.  (3)  from  a known  point-distribution 
along  the  surface  n=constant. 

Using  the  condition  for  local  orthogonality  (see  Appendix  A) , 

i.  e.  , 
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the  transformation  equations  (4-7)  can  be  written  as 

a «}.  t c + J^(Q  + B.  = 0 

nn  cc  n 5 

arid 

az  +cz__  + J^(Qz  +Rz)  = 0 

n T)  c £;  ^ 5 


(21) 


(22) 


,&t  the  surface  ri=constant,  all  derivatives  with  respect  to  c are 
easily  evaluated  from  the  given  point-distribution.  Therefore, 
if  the  n-derivatives  appearing  in  the  transformation  equations 
can  he  estimated  along  this  surface,  then  these  equations  can 
be  used  as  algebraic  equations  for  the  forcing  functions.  For 
a reasonably  fine  grid  in  the  n-direction,  as  would  normally  be 
used  in  a viscous  region  near  the  surface  n=constant,  it 
appears  reasoncible  to  assume  that 


z - 0 . (23) 

nn 

A good  estimate  of  and  is  obtained  by  linear  interpolation 

between  their  corresponding  values  at  the  end-boundaries  where 
these  quantities  are  known  from  the  specified  boundary-points. 
Finally,  z^  is  determined  from  the  orthogonality  condition  (20) 
as 

z = - /z  . (24) 

n n C C 

With  all  the  coordinate  derivatives  evaluated  along  the  surface, 

Q is  eliminated  between  Eqns.  (21)  and  (22)  to  yield  the 
following  expression  for  R; 


R = 


— [az  (J>  + c(z  (<>  - (ji  z )]  . (25) 

n nri  n CC  n cc 
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The  denominator  in  Eqn.  (25)  has  intentionally  not  been  repre~ 

sented  as  J^  in  order  to  admit  the  possibility  of  employing 

different  finite-difference  representations  for  the  derivatives 
2 

in  J and  in  (<|)  - <p^z  ) in  Eqn.  (25).  This  will  be  discussed 

n c c n 

further  in  the  next  section. 

Thus,  the  function  R can  be  determined  at  both  boundaries 
represented  by  Sind  n~n£  . Denoting  these  values  of  R by 

R^(i;)  and  R2(c)/  the  complete  function  R(n,i,)  to  be  vised  in 
the  interior  may  be  determined  as 

R(n,C)  = Rj^  + (R2  - f(n)  (26) 


where  f(n)  is  a normalized  monotonic  continuous  function  of  n» 
used  to  provide  a smooth  variation  for  R in  the  solution  domain. 
If  the  forcing  function  Q is  known  and  is  monotonic  with 
respect  to  H/  then  f(n)  may  be  chosen  as 

Q(n,5)  - Qi 

f(n)  = 1 . (27) 


Another  suitable  choice  for  f(n)  consists  of  the  sine  function 
or  the  exponential-sine  function. 

This  approach  has  been  used  to  determine  the  function  R 
from  the  specified  distribution  of  points  along  the  boundaries. 
The  resulting  solution  for  the  transformed  coordinates  confirms 
that  the  required  consistency  has  been  properly  maintained 
between  the  boundary  conditions  and  the  forcing  functions. 
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c.  Effect  of  Increasing  Magnitude  of  Forcing  Function  on 
Solution  Convergence  Rate 

It  has  been  reported  that  the  use  of  very  large  magnitudes 
for  the  forcing  functions  leads  to  numerical  difficulties  in 
terms  of  instabilities  and  oscillations.  However,  forcing 
functions  of  large  magnitudes  become  necessary  when  a highly 
non-uniform  coordinate  spacing  is  required,  such  as  in  a 
boundary  layer  for  a high~Re  flow.  A solution  to  this  situation 
frequently  consisted  of  obtaining  the  required  coordinates  by 
approaching  the  large  forcing  functions  in  a stepwise  increasing 
manner.  In  the  present  work,  it  became  possible  to  analyze 
the  cause  of  this  adverse  numerical  behavior  by  examining  the 
limiting  form  of  the  inverted  equations  at  the  boundaries. 

In  its  simplest  form,  this  limiting  equation  is  given  by  Eqn. 

(9)  for  a vertical  straight-line  boundary  represented  by 
;=constant.  Figure  5 shows  the  solution  of  Eqn.  (9)  for 
increasing  values  of  Q(n)  as  given  by  Eqn.  (12),  with 
equal  to  the  length  of  the  boundary  AB  in  Fig.  1.  Second-order 
accurate  central  finite-differences  have  been  used  to 
represent  all  derivatives  in  order  to  obtain  these  solutions. 

It  is  seen  from  Fig.  5 that  a satisfactory  solution  is 
obtained  with  l^lmax  — solution 

begins  to  violate  the  maximum  principle  and  shows  a clear 
oscillatory  behavior  for  =10^.  It  must  be  noted  that 

this  oscillatory  solution  is  truly  the  correct  converged  solution 
of  the  finite-difference  equation  used.  This  behavior  is  due, 
in  principle,  to  the  increasing  stiffness  of  Eqn.  (9)  with 
increase  in  Q.  The  situation  may  be  remedied  in  one  of  two 
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ways:^  the  use  of  a finer  mesh  in  the  transformed  coordinate  n 
or  the  use  of  directional  differencing.  The  use  of  a finer 
mesh  leads  to  an  increase  in  the  computer  storage  and  time 
required  for  the  solution  of  the  complete  coordinate  equations 
(4-7).  Therefore,  the  use  of  directional  o;;.  upwind  differences 
is  adopted  here  to  obtain  the  proper  solution,  with  large  O. 

For  this  purpose,  Eqn.  (9)  is  represented  in  its  linearized 
form  as 

^ + C,  - 0 (20) 

nn  -L  n 

and  = Q 4-^  • (29) 

The  term  41^  in  Eqn.  (29)  is  represented  by  a central-difference 
approximation,  while  41^  in  Eqn.  (28)  is  represented  by  an 
upwind  difference  such  that  a backward/forward  difference  is 
used  when  is  negative/positive.  This  is  easily  accomplished 
by  the  toilowing  generalized  representation: 


C,  4 - (l~u)  C,  4>, 


u 

5' 


ic,  4..  + c:  4..  i 


(30) 


'F 


■D 


f ^ 

if  central 

where 

U =“  c 

ll 

for  upwintl 

differencing  is  tu  be  used, 
dlf  fur  unci  nij  ; 


1 jcj 

and  4 , 4,,  ! 4,.  denote  tv/o-point  central-,  forward-  and 

C F U 

backv/ard"d.lf f crunce  roprusuntations , respectively, of  4^. 

Figure  S shows  tlie  effect  of  using  upwind  d if  lev  unci  ng  in  the 

solution  of  Eun.  (28).  With  IqI  in  the  vange  where  an 

' 'max 


acceptable  central-difference  solution  was  not  obtainable, 
upwind  differencing  leads  to  a completely  well-behaved  solution. 
Some  differences  are  also  observed  in  the  two  solutions  even 
for  |qL  ^ ■ 10^,  when  an  acceptable  central-difference  solution 
was  readily  obtainable.  Because  of  the  second-order  accuracy 
of  the  central-difference  solution,  versus  the  first-order 
accuracy  of  the  upwind-difference  solution,  the  former  may  be 
preferred.  However,  the  usefulness  of  upwind  differencing, 
with  large  Q,  indicates  a more  general  applicability  of  the 
upwind-difference  scheme. 

In  working  with  this  limiting  equation  (28) , it  was  also 
observed  that  the  central-difference  scheme  converged  very 
efficiently,  even  when  the  resulting  solution  was  not  monotonic. 
I'or  those  values  of  Q for  which  the  limiting  equation  (9)  leads 
to  a non-monotonic  solution,  the  solution  of  the  full  equations 
(4-7)  exhibits  instabilities  when  central  differences  are  used 
lor  all  the  derivatives.  On  the  other  hand,  upwind  differencing 
led  to  an  acceptable  solution  very  efficiently,  even  though  the 
convergence  crittirion  was  not  satisfied.  Efforts  were  made  to 
improve  the  convergence  properties  of  the  solution  by  employing 
ijuasi-linu.u  i/alion  in  Eqn.  (9),  i.e.,  by  representing  it  as 

* t (3Q  - (2Q  - 0 (31) 

wliere  the  quantifies  in  parentheses  are  evaluated  at  a previous 
iteration.  Uuasil inua- ion  was  found  to  increase  the  con- 
vergence rate  of  the  central-difference  scheme,  but  had  little; 
influence  on  the  upwind  sclieme.  The  following  table  uuimuari'/.es 
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the  effect  of  the  difference  scheme  on  the  convergence  behavior 


of  the  solution. 


Error  (e)  after  20 

Iterations 

of  Quasilinearized  Equation  (31) 

^Iraax 

Upwind 

Central 

10^ 

!-■ 

O 

1 

10-1^ 

lo" 

I-* 

o 

-14 

10  - oscillatory  solution 

ic6 

IMAX  n+i 

The  error  was  defined  as  e - I | (j) . where  n denotes 

i*=l  ^ 

iteration  number.  It  should  be  pointed  out  that,  even  with 

-4 

t “ 10  , the  upwind-difference  solution  was  a completely 

acceptable  solution  for  (|) . Therefore,  it  appears  that  appro- 
priate use  of  upwind  differencing  constitutes  an  efficient 
scheme  for  the  solution  of  the  complete  transformation  equations 
(4-7)  also.  In  the  interior  of  the  solution  domain,  it  is 
important  that,  if  upwind  differencing  j.s  used  for  it 
should  also  be  used  for  z ; similar  remarks  hold  for  (j)„  and  Zw. 
It  is  now  clear  that  upwind  differences  may  be  needed  to 
represent  the  first-order  derivatives  appearing  explicitly  in 

Eqn.  (25)  for  the  determination  of  R.  Due  to  the  different 

2 

origin  of  the  term  J'‘  in  Kqn.  (25)  , the  derivatives  appearing 
2 

in  J may  be  represented  by  central-difference  approximations. 
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d.  Determination  of  Forcing  Function  Q for  Boundary-Layer- 


Dependent  Mesh  System 

The  basic  numerical  transformation  procedure  provides  a 
means  for  generating  surface-oriented  coordinates,  with  the 
forcing  functions  providing  control  on  the  spacing  of  the 
coordinates  in  the  physical  domain.  However,  a systematic 
procedure  is  not  available  for  setting  up  the  variable  spacing 
of  the  coordinates  or  for  determining  the  forcing  function 
necessary  to  achieve  that  spacing.  In  the  present  work,  a method 
has  been  developed  to  compute  the  forcing  function  Q(n)  that  leads 
to  a mesh-point  distribution  along  a c-constant  coordinate  such  that 
certain  prescribed  criteria  are  satisfied.  The  boundary  layer 
theory  is  used  to  establish  an  important  one  of  these  criteria; 
hence,  the  resulting  mesh  system  is  referred  to  as  a boundary- 
layer  dependent  coordinate  system. 

The  discussion  presented  here  is  with  reference  to  flow 
past  an  airfoil.  However,  the  concept  is  general  and  can  be 
easily  adapted  to  other  flow  problems.  Based  on  experimentation 
with  an  inviscid  subsonic  flow  solution  for  an  airfoil  at 
incidence  of  up  to  30°  and  an  analysis  of  truncation  errors  of 
the  finite-difference  transformation  equations,  the  desired 
spacing  of  the  n-coordinates  is  required  to  satisfy  the  following 
criteria; 

a.  A minimum  number  K of  n-lines  must  be  placed  in  the 
boundary  layer  on  the  airfoil  surface . 

b.  The  maximum  step  size  (AiJ>)j,^ax  adjacent  to  the  far- 
field  boundary  should  be  approximately  unity  (in 
non-dimensional  variables)  except  when  this  boundary 
is  placed  at  infinity. 


26 


c.  The  grid  spacing  in  the  boundary  layer  should  corre- 
spond to  constant  increments  in  the  flow  variable 
with  the  largest  normal  gradient. 

d.  Second-order  derivatives  of  the  mesh  spacing  must 
be  minimized. 

The  first  two  criteria  serve  to  provide  appropriate  resolu- 
tion of  the  boundary- layer  region  near  the  airfoil  as  well  as 
the  far-field  region  of  the  flow  field.  Presently,  the  flow  in 
the  boundary-layer  region  is  being  represented  by  the  Blasius 
series,  for  the  purpose  of  the  mesh  generation.  Although  a more 
accurate  representation  of  this  flow  is  easily  obtainable,  the 
simple  Blasius  solution  was  used  in  developing  this  procedure. 
The  far-field  boundary  is  located  at  approximately  ten  chord- 
lengths  from  the  airfoil  surface  and  is  circular  in  shape; 
a more  appropriate  shape  for  the  outer  boundary  is  described 
later  in  this  report,  A total  of  44  points  is  used  along  the 
n-direction. 

The  last  two  criteria  are  used  to  minimize  the  truncation 
errors  in  the  finite-difference  solution.  Difference  approxi- 
mations for  the  Navier-Stokes  equations  are  usually  second-order 
accurate  in  the  transformed  coordinates.  Therefore,  the  trunca- 
tion error  in  these  approximations  is  proportional  to  third- 
and  higher-order  derivatives  of  the  flow  variables.  In  the 
boundary-layer  region,  normal  gradients  of  the  streamwise  velo- 
city w are  larger  than  most  other  gradients.  Therefore,  if 
the  consecutive  n-coordinates  are  placed  at  constant  increments 

4 

of  w,  as  suggested  by  Thompson  , then 

3w 

rr-  = constant 

■j  n 


and 


9^w  ^ 

9 n ^ 9ri  ^ 


= 0 
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Furthermore,  if  the  same  w-increraent  is  employed  for  the  n-lines 
at  each  streeunwise  location,  then 


2 3 

3W  3 W _ 3 W _ 

The  spacing  of  the  ^-coordinates  within  the  boundary-layer  can 
thus  be  determined  from  the  velocity  profiles  for  a boundary- 
layer  solution,  such  that  the  truncation  error  is  minimized. 

Since  the  analysis  essentially  optimizes  the  n-line  spacing 
with  respect  to  the  gradient  of  one  variable  in  one  direction, 
the  validity  of  the  analysis  depends  on  the  validity  of  the 
boundary- layer  approximation  for  the  problem  under  consideration. 

The  need  for  minimizing  the  second-order  derivatives  of  the 
coordinates  is  explained  as  follows.  For  a simple  stretching  trans- 
formation in  the  physical  plane,  the  truncation  error  associated  with 
a central-difference  approximation  for  a first-order  derivative  is 
obtained  as 

t<i!  ] w,,  + ... 

where  c denotes  central  difference.  Therefore,  the  truncation 

error  introduces  an  artificial-viscosity  term  proportional  to 

di  , so  that  (j)„„  must  be  minimized.  Moreover,  the  transformation 
n n n n 

of  the  second’-cr ier  derivative  as 


w.  “ [w  /<t)  ]„  “ 
9 n n c 


w. 


- a2 

- % 


w 


nn 


+ ^ 


nn 


w 


/ 


together  with  central-difference  representations  for  w^^  and 
shows  that  has  an  adverse  influence  on  the  diagonal  dominance 
of  the  finite-difference  equations.  This,  again,  indicates  that 
dij^,^  must  be  minimized. 
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All  of  the  above  criteria  must  be  satisfied  simultaneously. 
Also,  for  maximum  usefulness,  the  procedure  must  be  as  automated 
and  as  simple  as  possible.  Approximating  the  body  surface 
locally  as  a flat  plate,  the  transformation  equation  for  the 
normal  coordinate  is  given  by  Eqn.  (9) . With  a forcing  function 
of  the  form 

Q = -A  e"^'^  , (31) 

Eqn.  (9)  can  be  integrated  once  to  yield 

= 2 ^ e”^"^  . (32) 

d(p  0 

Here,  the  constant  of  integration  has  been  set  to  zero  by 
assuming  that  ari/34i  -►0  as  (ji  -►  41^.  Evaluating  Eqn.  (32)  at 
the  body  surface  n = 0 yields  an  estimate  for  A/D,  in  terms  of 
the  step  size  required  near  the  body  surface,  as 

§-^Un/iVn>'  • <”> 

Evaluating  Eqn.  (32)  at  ri  = (n^,„  " 4^)  yields  an  estimate 

luaX  ^ 

for  D,  in  terms  of  the  step  size  near  the  outer  boundary, 

as 

= '‘"12  I • (34) 

max 

To  place  K points  within  the  boundary .layer , the  forcing 


function  used  is 


where  nj^  ■ (k-l)An.  With  this  Q,  inteyration  of  Eqn.  (9)  leads 
to 


(ill) 


2 


where 


K 

I 


k=l 


-Dv{ri/n,  ) 
{e  ^ 


+ Isgn(nj^-n)  ,0]  [1-e  ^ ]} 

(36) 


M (n  »nj^) 


sgn(nj^“n) 


for 

n 1 \ 

• 

/ 

for 

n > nj. 

• 

t 

for 

Hk"^  1 

0 ; 

for 

nk-n  < 

0 . 

Evaluating  Eqn,  (36)  at  = (i+l/2)An  for  i = 1,  2,  (K-1)  , 
leads  to  (K-1)  simultaneous  equations  for  Aj^/D  for  k = 1,  2,  K 
with  D determined  from  a modification  of  Eqn.  (34)  as 


D 


= f s,n[2  g 


A /An 


A(|| 


2 

<)  ]/(n 


max 


max 


(37) 


where  f is  a weighting  factor,  - 1.1.  The  system  of  equations 
for  Aj^/D  is  completed  by  a constraint  relating  Eqns.  (32)  and 
(36)  so  that 
K 

I A,  = g A (38) 

k=l  ^ 

where  g is  a constant,  = 0.4. 

The  (K“l)  values  of  A<t'^  required  in  Eqn.  (36)  are  determined, 
so  as  to  minimize  9w/3n,  using  a Newton-Raphson  iteration  with 
the  Blasius  series  solution  for  the  velocity  E(n)  in  a flat-plate 
boundary  layer.  Here,  n is  the  boundary-layer  ordinate,  determined 
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— ■“'.~."'TOti’r"iimiT  iimii 


iteratively  by  the  relation 


(S.1) 


- (F^  - 2Wj^)/(F^  + 0.01) 


where  i = 1,  2,  K, 


F. 

1 


an,.  - 


a^  ~4 

rr  '’i 


lla‘ 

~Tr 


-7 

^i 


375a^  -10 


a = 1.32824  , F = 3F/8n  , = i/(K+l)  , 

and  s denotes  iteration  number. 

The  physical  ordinate  is  then  given  as 


(39) 


(fr.  „ = 2n.  /z/(w^ReT  (40) 

where  is  the  velocity  at  the  edge  of  the  boundary  layer.  If, 
instead  of  prescribing  K,  it  is  desired  to  specify  the  constant 
increment  Aw,  then  the  corresponding  value  of  K is  determined, 
by  integer  arithmetic,  using  the  relation 


K = (W  /Aw  - 0.01) . (41) 

e 

The  procedure  described  above  was  employed  to  generate  the 
forcing  function,  Q,  tor  various  values  of  the  parameters  for  a test 
case.  The  corresponding  (ji-distribution  is  then  determined  by  nu- 
merical solution  of  Eqn.  (9)  . The  second-order  derivative  was 
approximated  by  a central  difference,  while  upwind  differencing 
was  used  for  the  unknown  ({»^.  A quasilinearized  tridiagonal 
scheme  of  solution  requircsd  about  20  iterations , whereas  a 
linearized  successive  over-relaxation  (SOR)  scheme  with  optimized 
acceleration  parameters  required  about  40  iterations.  Figure  6 
shows  the  results  obtained  for  Q,  <{) , and  (|)^^  for  z = 1.0, 
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FORCING  FUNCTION  AND  COORDINATE  DERIVATIVES  FOR  BOUNDARV-LAYFP  DEPENDENT 
MESH  FOR  NACA0018  AIRFOIL,  Re  = 2 x 10®,  z = 1.0. 


Re  = 2 X 10^,  n = 43,  An  = 1/  and  W = 1.2,  Aw  = 0.15  so 
in^x  6 

that  K = 7.  The  computed  values  of  are  compared  with  the 
desired  values  of  in  Table  I which  also  lists  the  corresponding 
values  of  Q^.  It  is  seen  that  the  procedure  leads  to  a very 
satisfactory  mesh-point  distribution. 

The  procedure  has  also  been  applied  at  various  streamwise 
locations,  downs treeim  of  5 percent  chord,  on  a NACA  0018  airfoil. 
For  simplicity,  a Blasius  profile,  with  a constant  velocity 
at  the  edge  of  the  boundary  layer,  has  been  used.  The  first 
20  n-lines  in  a (71  x 44)  boundary- layer-dependent  mesh 
system  are  shown  in  Fig.  7 for  Re  = 2 x 10  , = 43,  An  = i# 

Wg  = 1.2,  AW  = 0.15,  with  K = 7. 

Some  improvements  can  be  made  in  this  automated  procedure , 
for  ex2unple/  by  using  the  appropriate  inviscid  velocity  and 
then  specifying  the  correct  boundary- layer  profile  from  a 
boundary- layer  solution.  However,  the  improvement  would  be 
limited  by  the  extent  of  validity  of  the  boundary-layer 
solution.  Therefore,  difficulty  may  be  encountered  in  regions 
of  large  surface-curvature  where  the  transformed  coordinate 
C=constant  may  not  be  orthogonal  to  the  body  surface. 
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TABLE  I 


SOLUTION 

OF 

(t>  4-  Q(j)^ 

NACA  0018 

= 0 WITH  BLASIUS  BOUNDARY- 

-LAYER  DEPENDENT 

Q (n) 

FOR 

AIRFOIL  AT  Re  = 2x10^,  z » 

1.0 

xl0“^ 

A . 

^i/Blasius 

/Computed 

0 

e. 72603 

0.000243 

0.000248 

1 

-0.19147 

0.000487 

0.000500 

2 

-0.27910 

0.000736 

0.000756 

3 

-0.67778 

0.000996 

0.001024 

4 

-1.1501 

0.001280 

0.001316 

5 

-1.7853 

0.001610 

0.001660 

6 

-1.8412 

0.002077 

0.002103 

7 

-1.1115 

0.002672 

8 

-0.67517 

0.003402 

9 

-0.40885 

0.004342 

10 

-0.24758 

0.005548 

11 

-0.14993 

0.007099 

12 

-0.09079 

0.009092 

13 

-0.05498 

0.011653 

14 

-0.03329 

0.014945 

15 

-0.02016 

0.019174 

16 

-0.01221 

0.024609 

17 

-0.00739 

0.031593 

18 

-0.00448 

0.040568 

19 

-0.00271 

0.052100 
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TABLE  I (contirued) 


xlO 


-6 


'ifOlaBiuB 


, computed 


i 

-0.00164 

•0.00099 

O.UUObO 

•0.00016 

O.OOOZZ 

•0.000,13 

-0.00000 

-0.0000,') 

■0.00003 

■0.0000?- 

■0.00001 

-0.00000 

-0.00000 

-0.00000 

-0.00000 

-0.00000 

-0.00000 

-0.00000 

-0.00000 

-0.00000 

-0.00000 

-0.00000 


0.066920 

o.oil*)'!!)! 

0,11043 

0.1410/ 

0. 

0 . 2 14 1 4 
0.  10l*/(> 
0.3tl62U 
0.4SI<)00 

0 , i/ibm 
0.016.')! 
1.0461 
,1.3361 

1. /O/'.’- 
2.1,696 
2.7416 
3.4370 
4.2631 
5.2183 
6.2904 
7.4596 
8.7033 


-0.00000 


10.0000 


1.  LOCATION  OF  FAF-FIFLD  UOUNDAUY  FOH  LimSONIC  IIXTFITNAL  FLOW 


a . Fac~Fiuld  Boundary  at  Ini'inity 

A major  quoutiun  that  ariuoa  in  thu  study  ot  subsonic  tlow 
past  a body  concurns  thu  nhapu  and  location  oi  the  iar-tiuld 
boundary  for  thu  solution  domain.  From  tho  viowpoint  of  analytical 
riqor,  this  Ixnindary  is  idoally  locatud  iniinituly  iar  irum  titu 
body,  Thu  unboundud  physical  domain  must  then  bu  mappud  to  a 
linitu  domain  by  an  appioprl  al »»  analytical  t ransiormation  bauod 
on  tliu  asymptotic  buhavior  ol  thu  Llow  lai  I rom  thu  body.  Thu 
analytical  tranr tormation  can  alrio  aui  vu  to  proviilu,  at  luant 
p.ait  ially,  thu  Mtiutcltiiuj  ol  soim*  ol  thn  critical  mnions  in  thu 
llow  iiuld.  Thu  ruBultinq  coordinatus  can  bu  iitadu  to  bu  surt'acu-' 
oriuntud  alonij  all  boundarius  uxcoiit  alonq  thn  surlacu  ol  thu 
arbitrary  body.  Thus^t  coordinatus  aiu  than  transtoriuud  nimiuri- 
oally  to  achieve  complutu  surlacu**oriuntatiun  and  thu  linal 
dusirud  strotchiny  oi.  thu  various  royions  ol  thu  problum  domain. 

Tltis  comi>ination  ol  analytical  and  mimorical  t.rnnslormatiuns  also 
considurably  ruliuvus  thu  burdun  oi  thu  numurical  transtormatioii . 

As  royards  thu  uhapu  ol  thu  lar-liuld  boundary,  it  appears 
that,  lor  a laiyo  class  ol  two-dlmunuional  or  axiuyiKunutric  Hows 
past  linitu  bodies,  a parabolic  shape  may  bu  suitable  lor  an 
outer  boundary.  This  coniiyurution  purmits  ono  ol  thu  translormud 
coordinates  to  bu  oriuntud  alony  tin?  same  yunural  diruotion  as 
the  main  ilow,  so  that  thu  associatud  numurical  solution  can  have 
bettur  convurguncu  propurtius.  Thuruloru,  it  may  bu  moru  appropriati 
to  uut  up  thu  numurical  translonuation  in  turms  ol  parabolic 
coordinatus  rathur  than  Cartesian  or  cylindrical  coordinates. 
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after  mapping  thorn  to  a finite  domain.  Thus,  tho  total  trann- 
formation  procedure  would  conaiet  of  tho  following  etepat 


I’hyuical 


coordinate 


ayutem:  {^,7.) 

I 


(li' 


i 


(N,t;) 


parabolic  trana formation 
(analytical) 


finite-region  mapping 
.analytical  contraction. 


I'inal  tranaformed 
coordinate  ayatuiu: 


I numerical  trana formation 

(iwc) 


'LVo  approachua  proaent  thumaulveu  at  thia  utage,  regarding  the 
ogu.il  iona  to  bo  luuui  to  accompliah  tho  numerii^al  trana for»iuition. 
in  the  eguationa  uaed  ao  far,  the  elliptic  operator  employed  haa 
been  tlie  baplacian  operator.  In  order  to  conform  to  thia,  it  would 
now  be  neceaaary  to  formulate  tho  baplacian  operator  in  termu  of 
the  (Njij)  coordinutea  aa  independent  variableu.  Tho  corroapondi ng 
I'oiaaon  oguationa  Lor  tho  coordinatou  n,  t,  then  nood  to  ho  invortoil 
to  make  t| , t.  tho  independent  variablou  and  N,  b tho  dependent 
variabloa.  ThoroLoro,  in  thia  approach,  tho  inuaorioal  tranufoniia- 
tion  would  employ  the  true  baplacian,  no  that  the  oquat  ionn  to  )m* 
inverted  are 


V 


2 

N,S 


<1 


0 


'N 


3y 


and 


- U 


(42a, b) 


wharo  ia  the  Laplacian  in  (N,S)  coordinates  and  is 

obtained  as 

..2 


N,S 


1 ' 5 ' ' n ' 9 ' ' a 

— y [ (N  ) — — K-  4-  N -Tm  ) “"5"  S ^q]  (43) 


with  primus  denoting  differentiation  of  a variable  with  respect  to 
its  argument.  Also,  h is  the  scale  factor  for  the  parabolic 
trans  formation 


“ + ^+p> 


(44) 


b'o  that 
h^ 


^2  + .2 
V V 


(45) 


The  inversion  of  Uqs,  (42a,b)  leads  to  the  following  equations  for 
the  numerically  transformed  coordinates: 


a N 4 

i|il 

2 b 

N , 4-  c N,., 
nc,  fA 

• 1 

“ N 

» 

4-  h^ 

j2 

4-  KN  ) 

. 0 

a S 4 

nil 

2£ 

S 4 c S .. 

DC.  (.C 

1 1 

. - s 

I 

4.  h2 

j2 

+ UK^) 

- 0 

(46a, b) 


where 


h - 


2 ' 9 2 '2 

U 4 N ^ 

-{N  N y*2  y s n'2) 

n c n C 


2 '2  2 '2 
N b 4-  U N 

II  n 


NS-  NS 
I)  C CM 


(47) 
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2 (III 

It  should  be  noted  that  h as  well  as  the  quantities  N , N , S , 

% I 

S must  be  expressed  in  terms  of  N and  S in  the  above  equations. 

The  details  of  the  derivation  of  Eqs,  (46a,b)  are  described  in 
Appendix  B. 

In  an  alternate  approach,  it  seems  that,  for  the 
coordinate  equations,  perhaps  it  is  not  essential  to  work  with  a 
true  Laplacian  operator.  Another  elliptic  operator  that  satisfies 
the  maximum  principle  would  also  bo  satisfactory.  Of  course, 
the  choice  of  the  Laplacian  was  guided  by  its  similarity  to  the 
equations  for  the  corresponding  inviscid-f low  stream  function 
and  potential  function.  For  reasons  of  mathematical  simplicity, 
the  following  system  of  equations  for  the  numerically  transformed 
coordinates  may  be  preferable: 

*'nN  ’^SS  " ^ ’ ‘’NN  *'SS  " ^ • (48a, b) 

The  corresponding  inverted  equations  for  the  N,S  coordinates  are 
of  the  same  fom\  as  Eqs.  (4-7). 

Boundary  conditions  specified  in  the  ((j),z)  domain  must  be 
expressed  in  terms  of  the  (N,S)  variables,  through  use  of  the 
parabolic  transformation,  followed  by  the  analytical  contraction 
trans  formation . 

Before  describing  any  results  obtained  with  these  combined 
transformation  procedures,  the  discussion  presented  in  the  following 
section  must  be  taken  into  consideration. 
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b. 


Far-Field  Boxindary  at  Finite  Distance  From  Body 

In  many  practical  situations/  it  may  not  be  essential  to 

locate  the  far-field  boxandary  at  true  infinity.  For  example, 

analyses  of  flow  over  a cylinder  have  frequently  employed  a far- 

field  boundary  located  at  2d)out  10  cylinder-radii  from  the 
I 9 

cylinder  . As  shown  below,  this  is  the  location  where  the  free- 
stream  pressure  is  disturbed  by  less  than  one  percent  due  to 
the  presence  of  the  cylinder  in  the  flow.  For  incompressible 
inviscid  flow  past  a cylinder  of  radius  r^,  the  pressure  coef- 
ficient Cp  is  expressed  as 


where  the  subscript  « denotes  f ree-strecim  conditions  . Here , V is 
the  local  total  velocity  in  the  flow  given  as  follows,  in  terms  of  the 
radial  and  tangential  velocity  components  u^  and  u^,  respectively, 


V 


2 


+ u 


2 

0 


where  , 

u^  = u^(i  - -|)  cose 

J.2 

^0  “ sine 


(50) 


(51) 


Substituting  Eqs.  (50-51)  in  Eq. 

P " 2 (r/r^)^ 

This  implies  that 


(49)  gives,  after  simplification, 

(52) 


max 


r/r. 


10  . 


0.01 


when 


(53) 


Therefore,  an  appropriate  far-field  boundary  for  the  flow  past 
a cylinder  would  encompass  at  least  the  region  shown  in  the  fol- 
lowing sketch. 


FIGi  8,  LOCATION  OF  FAR-FIELD  BOUNDARY  FOR  FLOW  PAST  A CIRCULAR 
CYLINDER. 

For  viscous  flow  calculations,  it  would  be  more  desirable  to  have 
a far-field  boundary  such  that  the  distance  A'C  increases 
approximately  as  the  square-root  of  the  streamwise  distance. 

The  above  analysis  can  also  be  employed  to  estimate  the  extent 
of  the  region  in  which  the  inviscid  pressure  for  flow  past  an 
airfoil  is  disturbed  by  less  than  one  percent.  A symmetric 
Joukowski  airfoil  is  considered  here  for  this  purpose.  With  a 
Joukowski  transformation,  a cylinder  of  radius  (r^  + a) , with 
center  at  (-a,  0) , transforms  to  a Joukowski  airfoil  of  chord 

n 

SL  s 4r^(l+E  ) and  thickness  ratio  t - (3/3’/4)e,  where  e is 

the  eccentricity  defined  as  e = Since  the  inviscid  flow 

past  the  airfoil  is  obtainable  by  transformation  of  the  flow  past 
a cylinder,  the  region  of  disturbed  inviscid  flow  for  the  airfoil 
is  also  obtainable  by  transformation  of  the  region  of  disturbed 
inviscid  flow  for  the  cylinder.  Accordingly,  it  is  necessary  to 
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transConn  a larger  circle  of  radius  (lOr^  + 10a)  with  center  at 

(-10 a,  0)  using  the  Joukowski  transformation.  This  results  in 

2 

a larger  Joukowski  airfoil  of  chord  L “ 40  r (1  + e ) and 

o 

maximum  thickness  “ 30/3"  a.  The  eccentricity  and  the 

max  ■* 

thickness  ratio  of  the  larger  airfoil  are  the  same  as  those  for 
the  original  aivfoil.  With  reference  to  the  chord  of  the  original 
airfoil,  the  ^joundary  of  the  region  of  disturbed  inviscid  flow 
past  the  airfoil  is  obtained  as  a large  Joukowski  airfoil  of 

chord  L/i,  = 10  , 

15  nr 

maximxam  thickness  T,,„/A  “ /3  e 13e  . (54) 

IliaX  ^ 


Therefore,  for  a 10-percent  thick  airfoil  at  zero  incidence,  a 
far-field  boundary  based  on  the  region  of  1-percent  inviscid- 
pressure  disturbance  would  appear  as  shown  in  the  following 
sketch. 
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by  (10  sine)  to  account  for  the  increase  in  effective  frontal  area 
of  the  airfoil. 

It  is  necessary  to  observe  that  the  above  analysis  only 
provides  the  location  of  the  contour  of  1 Cp  1 = 0.01. 

This  serves  as  a useful  yuideline  for  the  location  of  the  far- 
field  boundairy.  However,  by  employing  the  undisturbed  free-otream 
conditions  at  the  location  of  1-percent  disturbance,  the  error 
introduced  in  the  results  at  the  body  surface  may  exceed  1 percent 

A study  was  made  to  assess  the  inferences  of  this  analysis 
by  computing  the  inviscid  flow  past  a symmetric  Joukowski  airfoil 
of  approximately  9-percent  thickness  ratio  (i.e.,  t „/s,  = 0.09), 
at  various  angles  of  incidence.  liie  outer  boundary  for  these 
calculations  was  an  ellipse  of  major  axis  2z_^  and  minor  axis  2^^, 
with  the  major  axis  aligned  with  the  free-stream  direction.  The 
results  of  these  calculations  are  summarized  in  Table  II . Compari 
son  is  made  of  the  coefficients  of  minimum  pressure  , 

drag  lift  and  moment  about  the  leading  edge 

various  values  of  z and  (b  . For  each  value  of  S studied,  the 
last  row  corresponds  to  the  values  of  z^  and  (|)^  obtained  on 
the  basis  of  the  above  analysis.  It  is  seen  that,  with  the  far- 
field  boundary  located  in  accordance  with  the  prediction  of  the 
foregoing  analysis,  the  maximvim  error  in  and 

1 percent.  The  maximum  error  in  is  about  0.02.  However,  the 

maximum  error  in  C , . s is  almost  2.7  percent  and  occurs  for 

p (min) 

B = 10®  and  30°.  But  it  should  be  noted  that,  when  has 

this  large  error  for  the  smallest  z^,  (f)^  investigated,  a major 
part  of  this  error  has  been  already  incurred  at  the  largest 

(j)  investigated. 
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TABLE  II  I 

EFFECT  OF  OUTER-BOUNDARY  LOCATION  ON  SURFACE  COEFFICIENTS 

FOR  INVISCID  FLOW  PAST  A SYMMETRIC  JOUKOWSKI  AIRFOIL  | 
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0 

2 /C 

Op' 

'“p  ,min 

‘^M,LE 

1 

• 

OP 

OP 

-0.34513 

0.0 

0.0 

o 

• 

o 

] 

20 

20 

-0.34610 

-0.205  xlO"^ 

-0.820  XIO”^ 

-0.05  xlO"^ 

1 

0“ 

10 

10 

-0 .34616 

-0.196  xlO"^ 

-0.083  xl0“^ 

+0.21  xlO"^ 

1 

I 

10 

6 

-0.34640 

-0.189  xlO"^ 

-0.069  xlo"^ 

+0.16  xlO”^ 

1 

10 

4 

-0.34700 

-0.179  xl0“^ 

-0.290  xlo”^ 

-0.05  xlO”^ 

t 

5 

2 

-0.35060 

-0.159  xlO"^ 

-0  .163  xlO"^ 

-0.22  xlO"^ 

1 

• 

\ 

«e 

OP 

-9.2250 

0.0 

1.1667 

«0.2917 

1 . 

( 

20 

20 

-9.4196 

-0.0028 

1.1700 

0 .2922 

' 

1 

10® 

10 

10 

-9 .4190 

-0.0027 

1.1708 

0 .2924 

1 

1 

10 

6 

-9.4200 

-0.0024 

1.1727 

0.2932 

1 

10 

4 

-9.4390 

-0.0022 

1.1770 

0 .2949 

5 

4 

-9.4750 

-0.0029 

1.1780 

0 . 2948 

\ 

PP 

PP 

-75.550 

1 

o 

• 

o 

3.3593 

•=<0.8398 

1 

1 

30® 

20 

20 

-77.606 

-0.0186 

3.3790 

0.8448 

10 

10 

-77.594 

-0.0187 

3.3800 

0.8451 

This  indicates  that,  at  least  from  inviscid-flow  considera- 
tions, the  far-field  boundary  may  enclose  a considerably  smaller 
region  than  employee  so  far  (z^  = = 10  for  0=  0),  without 

significantly  affecting  the  accuracy  of  the  results.  The  con- 
clusion arrived  at  in  Ref.  9 for  the  proper  values  of  z , <J) 
must  be  because  of  the  considerably  larger  values  computed 

■k 

there  for  AC,,  and  AC^  for  0 = 0.  With  the  general  procedure 
very  similar  to  that  in  Ref.  9,  the  present  results  shov/  much 
smaller  AC^  and  AC- . A closer  examination  reveals  a difference 

iJ  1j 

in  the  manner  C is  integrated  along  the  body  surface,  to  compute 
P 

Cr^  and  Ct  • Considering  the  trapezoidal  rule,  the  integral  of 
D u 

Cp  around  the  body  contour  was  represented  in  Ref.  9 as 


IMAX 


# c ds  = Cf  c g as  - I [f.  . f..,i 


i=2 


(55) 


where  f.  = C g.  and  g consists  of  the  scale  facotrs 

Pi  1 

of  the  transformation  evaluated  at  the  body  surface  s. 

However,  since  g enters  this  integre  through  the  relation 
ds  = g d?  , it  has  been  found  more  appropriate  to  represent  the 
integral  in  the  following  form; 


^ C ds 
P 


IMAX 


I 


i=2 


I 


(56) 


Experience  with  other  problems^^  has  also  indicated  that  the 
product  of  terms  in  integrands  be  represented  as  the  product  of 
the  average  values  of  the  terms,  rather  than  as  the  average 
of  the  product. 

* Here,  the  prefix  A denotes  deviation  from  the  theoretical  values. 
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3.  ORTHOGONALIZATION  OF  COORDINATES 

As  mentioned  in  the  Introduction,  it  is  frequently  desirable 

to  work  with  an  orthogonal  or  near-orthogonal  coordinate  system. 

Non-orthogonality  leads  to  a non-zero  coefficient  of  the  mixed- 

derivative  terms  in  the  coordinate  equations  as  well  as  in  the 

flow  equations.  Most  numerical  schemes  available,  even  implicit 

ones,  treat  mixed  derivatives  in  an  explicit  manner.  This  tends 

to  retard  the  convergence  rate  of  the  numerical  solution. 

Mitchell^^  has  discussed  a second-order  accurate  finite-difference 

representation  for  treatment  of  the  mixed  derivative.  The 

representation  takes  different  forms  depending  on  the  relative 

values  of  the  coefficients  a,  b,  c in  the  governing  equations 

12 

(4-5) , The  further  analysis  of  Greenspan  and  Jain  relaxed 
some  of  the  conditions  obtained  by  Mitchell.  However,  implementing 
these  representations  in  the  implicit  alternating-direction 
numerical  scheme  presently  emp'loyed  did  not  lead  to  a gain  in 
computational  efficiency.  It  is  possible  that  this  implementa- 
tion of  the  mixed  derivative  was  still  not  fully  implicit; 
further  work  is  indicated  in  this  area.  These  considerations 
do  not  arise  in  the  case  of  nearly  orthogonal  coordinates,  Wiien 
the  mixed-derivative  terms  have  coefficients  that  are  nearly 
zero. 

An  effort  is  made  in  the  present  work  to  develop  a pro- 
cedure for  obtaining  a system  of  orthogonal  or  near-orthogonal 
coordinates,  starting  with  a given  non-orthogonal  coordinate 
system.  The  starting  system  is  determined  as  the  solution  of 
the  coordinate  equations  (4-7)  using  the  appropriate  forcing 
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I 
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function  Q for  the  end  boundaries,  with  consistent  distribution 
of  the  points  along  these  boundaries.  The  forcing  function  R 
must  be  consistent  with  the  desired  distribution  of  points  on 
the  lower  and  the  upper  boundaries.  Also,  the  boixndary 
coordinates  must  be  r-rthogonal  at  the  corners. 

This  solution  generally  yields  a good  set  of  body-oriented 
coordinates.  However,  the  coordinate-system  v-rill  not  necessarily 
be  orthogonal  because  the  distribution  of  points  on  the  far- 
field  boundary  may  not  be  consistent,  as  regard  orthogonality, 
with  the  point  distribution  on  the  lower  boundary  containing 
the  body  surface.  Therefore,  the  coordinates  must  be  adjusted 
so  as  to  produce  an  orthogonal  system.  For  this  purpose,  the 
streamwise  coordinates,  i.e.,  the  curves  of  n^constauit , are 
retained.  The  c-coordinates  are  altered  so  as  to  satisfy  the 
condition  of  local  orthogonality,  Eq.  (20),  i.e. 

+ 2^  = 0 . (20) 

Starting  at  the  lower  boundary,  the  c-curves  are  altered  in  a 
step-by-step  manner,  over  each  An  step,  as  illustrated  below. 


FIG.  10.  STEP-BY'STEP  ORTHOGONAL I ZAT I ON . 


.*1 
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The  above  sketch  shows  two  coordinate  curves  n ==  n-j^  and 
n = + Ati/  with  dots  denoting  the  starting  position  of  the 

points  obtained  from  the  solution  of  the  coordinate  equations . 
It  is  reqxiired  to  determine  the  positions  of  points  such  as 

I I I 

b , c , d , etc.,  such  that  the  resulting  coordinate  system  is 
orthogonal.  To  locate  b , for  example,  (}>^  and  atre  computed 
as  the  average  of  their  values  at  points  '2'  and  b,  using  the 
non-orthogonal  solution  along  (n^^+ATi)  . Thus, 


2b 


1 ,^'*>3“‘*’l^■^^^c"‘^a^ 

1 ^ 5ac ^ 


and 


2b 


1 ,<^3-^l^  + ^^c-^a^ 

2 ^ TKl ^ 


(57) 


Next,  and  z^  are  computed  from  the  known  non-ortho gonal 
solution  as 


2b 


At! 


V"2 

An 


(58) 


Then,  the  left-hand  side  of  the  condition  of  orthogonality, 
Eq.  (20)  , is  evaluated  to  yield  62jj  as 


z 

n 


2b 


(59) 


If  621^  ^ 0,  the  calculation  proceeds  to  the  next  point  along  the 
(n^^+An) -curve,  to  evaluate  520'  ^2d'  until  a change  in  sign 

is  encountered  by  the  left-hand  side  of  Eq.  (20) . For  the  sketch 

I 

shown,  this  would  occur  between  b and  c for  point  '2*.  Then,  b 
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is  located  between  b and  c,  using  linear  interpolation,  yielding 
62Jji  = 0 and,  hence,  a locally  orthogonal  coordinate  system 
in  the  vicinity  of  point  '2'. 

The  procedure  is  repeated  for  point  '3'  on  the  n^^-curve, 
beginning  with  point  b on  the  (rij^+An)  curve  and,  similarly,  for 
all  points  up  to  ) on  the  n^^-curve.  The  procedure  is 

then  repeated  for  the  pair  of  n-coordinates  given  by  (nj^+An)  = 

I I I 

constant  and  (nj^+2An) “constant,  using  the  points  b , c , d , 
etc.,  on  the  (n^+An) -curve.  Continuing  successively  for  all 
the  Ti-coordinates  until  the  upper  boundary  ’^“’^inax 
encountered,  this  will  result  in  a coordinate  system  that  is 
locally  orthogonal  at  all  the  grid  points.  It  is  to  be 
noted  that  only  two-dimensional  configurations  can  be  considered 
by  this  approach,  at  present. 

The  procedure  performs  best  when  the  n -coordinates  have 
only  convex  curvature  or  minimal  concave  curvature  and,  in 
these  cases , o restriction  is  required  on  the  degree  of  non- 
orthogonality of  the  starting  coordinate  system.  However,  when 
the  n-coordinates  contain  regions  of  large  concave  curvature, 
the  procedure  may  lead  to  cross-over  of  the  new  c-coordinates 
unless  the  local  mesh  aspect  ratio  (Az/A(}>)  in  the  physical  plane 
equals  or  exceeds  unity.  This  requirement  may  be  related  to 
the  Cauchy-Riemann  conditions  which  are  more  restrictive  than 
the  condition  of  orthogonality  [see  Appendix  A] . 
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a.  Partial  Orthogonalization 

Situations  may  arise  where  the  surface  coordinate  n " 
may  contain  a point  of  discontinuous  slope,  e.g.,  near  the 
trailing  edge  of  an  airfoil.  For  these  configurations,  it  is 
desirable  to  require  the  altered  coordinate  system  to  be 
orthogonal  everywhere  except  in  the  vicinity  of  the  point  of 
discontinuous  slope.  In  this  vicinity,  the  coordinate  curves 
are  required  to  intersect  at  a specified  angle  related  to  the 
angle  in  the  surface  coordinate  at  this  point.  This  may  bo 
achieved  by  a procedure  described  below. 

The  following  sketch  shows  a portion  of  the  surface 
coordinate  n ■ which  has  a discontinuity  at  point  C.  This 
discontinuity  may  be  represented  by  the  angle  y,  which  differs 
from  IT . 


In  terms  of  surface  points  immediately  adjacent  to  point  C, 

Y = TT  - tan“^  \b,^/Lz\  . ^60) 

For  partial  orthogonalization  near  the  point  C,  the  angle  between 
the  n-coordinate  line  through  C and  the  n^^-coordinate  should 
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be  y/2  on  either  side  of  the  discontinuity.  Now,  the  angle 
of  intersection  between  the  n“  and  ^-coordinates  can  be  obtained 
from  the  relation 


tana  * tan (a  -a  ) 
n C 


tana  - tana 

^ ^0 s;— 

1 + tana  tano 
n C 


(61) 


Noting  that  tana^ 
Eq.  (60)  becomes 


li-l 

DZ  I C 


(b  /'L  and  tana  “ = A /z  , 

n "n  32  n C 


tana 


i.e. 


tanc 


d)  2 - A Z 

C 

2 2 + di  4 

c n ^c^ti 


(62) 


Using  Eqs.  (6-7) , this  can  be  expressed  simply  as 

tana  “ J/b  . (63) 

Since  b is  a measure  of  non-orthogonality  of  the  coordinates, 

Eq.  (63)  implies  that  b approaches  2ero  where  the  coordinates  are 
nearly  orthogonal.  On  the  other  hand,  the  Jacobian  J is  never 
2ero  or  infinite  for  a one-to-one  transformation.  Therefore,  it 
is  better  to  work  with  cot a expressed  from  Eq.  (63)  as 

cota  = b/J  . (6  4) 


Therefore,  the  partial-orthogonalization  procedure  consists  of 
defining  a^  j,  for  each  computational  point  (i,j)  as 
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f : 


cota^^j  ■ - [8gn(  j-jc)](cot  e 


-f.  . -ki^j 


■i/ j 


cos  f.  . cos  g.  . 
1 / j 1 » j 

(65) 


where 


sgn (j-jc) 


( -1  for  (j-jc)  ^ 0 , 

\ +1  for  (j-jc)  > 0 , 


^ ‘‘’lmax,j  " n,jc 


’iij 


1 [ ’^j'/j  “ 

2 2.  J_  - z. 


'i/ jo 


‘1/ jc 


with 


{‘ 

^ Jmax 


for 

for 


j < jo 
j > jc 


(66) 


The  cosine  functions  have  been  included  in  Eq.  (65)  in  order  to 
ensure  complete  orthogonality  near  the  remaining  three  boundaries 
of  the  problem.  Use  of  the  function  [sgn( j-jc) ] ensures  that  the 
curves  of  ^“constant  are  nearly  parallel  in  the  vicinity  of  the 
point  G.  Here,  jc  represents  the  value  of  j at  the  point  c. 

Hereafter,  the  procedure  is  similar  to  that  described  in 
the  preceding  section  for  complete  orthogonalization , with  Eq. 

( 59)  now  being  replaced  by  the  following  for  partial  orthogonali- 
zation: 


'P,2b  * - (b/J)2b 


(67) 


where  6p  2]^  is  the  error  in  Eq.  (64)  if  the  C -coordinate  at  point 
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j'l 


'2'  also  passes  through  point  b.  In  fact,  Eq.  (59)  may  be  viewed 

as  a special  case  of  Eq.  (67)  when  cota.  . =0,  so  that  the 

1 # J 

complete-orthogonalization  procedure  is  a specialization  of  the 

partial-orthogonalization  procedure . 

Finally,  it  should  be  recognized  that  the  complete- 

orthogonalization  procedure  outlined  in  the  preceding  section 

13 

is  equivalent  to  that  of  Potter  and  Tuttle  who  achieve 

orthogonalization  by  solving  a Neumann  boundary-value  problem 

between  each  two  pairs  of  n=constant  curves.  It  should  be 

possible  to  generalize  the  work  of  Ref.  13  to  include  partial- 

orthogonalization  by  modifying  the  homogeneous  Neumann  boundary 

conditions  to  nonhomogeneous  Neumann  boundary  conditions  involving 

cota,  . as  defined  in  Eq.  (65)* 

1 / 3 

SECTION  IV 

OPTIMIZATION  OF  NUMERICAL  METHOD  OF  SOLUTION 

The  discussion  in  all  of  section  3 was  directed  towards 
obtaining  a suitable  coordinate  system,  using  as  much  formalism 
in  the  procedure  as  possible,  so  as  not  to  unduly  increase  the 
overall  computer  time  required  for  generating  the  desired  mesh. 

For  example,  the  use  of  upwind  differencing  in  generating  highly 
non-uniform  mesh  spacing,  as  well  as  the  consideration  of  an 
analytical  coordinate  transformation  simultaneously  with  the 
numerical  transformation,  may  be  also  viewed  as  optimization  of 
the  numerical  method  of  solution  of  the  transformation  equations. 
In  addition,  however,  some  measures  have  also  been  taken  to 


J 
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directly  optimize  the  numerical  method  itself.  Two  such  measures 
have  been  used  previously  and  are  described  elsewhere  . These 
consist  of  the  use  of  an  implicit  numerical  scheme  and  the 
simultaneous  solution  of  the  coupled  equations  (4)  and  (5)  for 
the  transformation.  Another  important  factor  contributing  to 
the  efficiency  of  the  nimierical  solution  is  the  starting 
solution  used  for  these  nonlinear  equations;  this  is  described 
here . 

1.  INITIALIZATION  PROCEDURE 

Frequently,  no  special  care  is  taken  to  formulate  the 

starting  conditions  used  in  an  iterative  numerical  solution 

procedure.  This  is  based  on  the  argument  that  initial  errors 

decay  with  exponential  rapidity  and  have  no  significant 

influence  on  the  solution  convergence.  However,  for  complex 

problems  with  severe  nonlinearities,  this  has  been  often  found 

to  be  untrue.  In  these  situations,  true  nonlinear  instability 

has  been  observed,  i.e.,  small  initial  errors  decay  whereas 

large  initial  errors  cimplify  and  cause  the  solution  to  become 

unbounded.  Therefore,  in  these  cases,  the  starting  solution 

employed  plays  a significant  role  in  the  success  of  the  overall 
14 

solution  . In  general,  the  solution  of  some  simplified 

limiting  form  of  the  governing  differential  equations  provides 

15 

a good  starting  solution  . However,  for  the  coordinate 
equations,  even  this  may  not  be  necessary.  The  following  simple 
geometrical  considerations  lead  to  a suitable  initial  distri- 
bution for  the  coordinates. 
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If  the  initialization  proceeds  in  the  direction  of  ih~ 
creasing  counters  i and  j , the  increment  along  a coordinate  at 
a given  station  can  be  made  to  have  the  same  ratio  to  the  total 
interval  along  that  coordinate  at  that  station  as  the  corresponding 
ratio  at  the  preceding  station.  This  permits  explicit  calculation 
of  the  initial  condition.  Thus,  with  reference  to  the  sketch 
shown  in  Fig.  12, 


JljJl 


- *1-1 . 


‘*'lMAX,j  "'*’l,j 


'^i,j-l  ~ ‘‘*1-1, j-1^ 
*IMAX,j-l  "***!, j-1 


z . . - z . . , 

-■Lfl  ■ 

''i,JMAX  ”^i,l 


“i-l,JMAX  "^i-1,1 
(68a, b) 


Similarly,  it  is  possible  to  formulate  analogous  expressions  for 
initial  values  based  entirely  on  information  at  a station  following 
the  station  under  consideration  and  set  up  the  initialization  to 
proceed  in  the  direction  of  decreasing  counters  i and  j.  Accordingly, 


**’i,j  ~ **’itl,j  ^ '**i,j+l  ~ *^1+1, i+1^  ^i,-]  ~ ^i,i+l  _ ^i-H,j  ~ ^i+l,jtl 

'^IMAX,j  "**’l,j  '*'lMAX,j+l  “‘*’l,j+l'  ^i,JMAX  “’^i,!  '^i+l,JMAX  '’“i+1,1 

(69a, b) 


Ideally,  the  initial  values  used  should  be  based  on  informa- 
tion at  stations  preceding  as  well  as  following  the  station  under 

consideration.  Thus,  (J> . . should  be  initialized  using  a suitable 

1 / J 

weighted  combination  of  relations  (68a,  69a)  and  z.  . using 

^ f J 

relations  (68b,  69b) . However,  only  one  of  the  right-hand 
quantities  would  be  known  in  this  combination,  depending  on  the 
direction  in  which  the  initialization  proceeds  with  respect  to 
i and  j.  To  maintain  the  initialization  procedure  non-iterative, 
the  unknown  quantity  must  be  approximated  in  a suitable  manner. 
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This  leads  to  the  following  simple  expressions , written  here  for 
increasing  i and  j. 


~ ‘^i-l,j 
" ‘^l,j 


and 

z . . - z . . , 

_ i>3-l 

^i,J.MAX  "^i,l 


^ “^IMAX^j-l  ^l,j-l  ^ 


’*’i,JMAX  ~ ‘^i-l,JMAXj 
'*’lMAX,JMAX  "‘*’l,JMAX 


(70) 


i-l,JMAX  i-1,1  IMAX,7MAX  IMAX,1 

(71) 


where 


w = ^i-l.JMAX  ~ ^i-l,j  . ^ ^ ^IMAX>j-l  ~ =^i,j-l  ^ 

‘^i-l,JMAX  " **’i-l,j-l  ' ^ ^IMAX,j-l  " ^i-l,j-l 

(72) 


This  procedure  has  been  employed  to  provide  the  starting 
solution  for  a variety  of  configurations,  in  order  to  test  it  for 
versatility.  For  a configuration  involving  a simple  rectangular 
domain  but  requiring  non-uniform  grid-point  distribution,  Eqs . 
(70-72)  provide  the  final  solution  exactly,  consistent  with  the 
grid-point  distribution  specified  at  the  boundaries.  Therefore, 
for  this  configuration,  the  converged  solution  of  the  transfor- 
mation equations  is  obtained  in  one  iteration  (Fia.  13). 

For  a typical  turbine-cascade  configuration,  the  starting 
solution  provided  by  Eqs.  (70-72)  is  shown  in  Fig.  14a;  the 
corresponding  converged  solution  is  shown  in  Fig.  14b.  The 
striking  resemblance  between  the  starting  solution  and  the  final 
solution  is  responsible  for  the  rapid  convergence  of  the  numerical 
method.  In  fact,  this  pai'ticular  configuration  was  used  as  a 
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FIG.  12.  INITIALIZATION  PROCEDURE 


ITERATION  NO. 


INITIAL  CONDITIONS  AND  FINAL  SOLUTION  FOR  COORDINATES  FOR  FLOW  THROUGH 
A STRAIGHT  CHANNEL. 


test  case  to  confirm  the  nonlinear-instability  phenomenon  caused 
by  large  initial  errors.  With  starting  conditions  consisting  of 
constant  values  of  (}(  and  z in  the  interior,  the  solution  failed 
to  converge. 

The  usefulr 2ss  of  Eqs.  {70-72)  for  an  external  flow  con- 
figuration is  indicated  in  Fig.  15  which  shows  the  starting 
solution  for  the  coordinates  suggested  for  analyzing  the  symmetric 
flow  past  an  airfoil.  Other  configurations  have  also  been 

g 

tested  by  Plant  . 

It  is  important  to  note  that,  in  certain  limiting  situations 
involving  vertical  or  horizontal  straight-line  boundaries,  care 
must  be  taken  to  avoid  possible  division  by  zero  in  Eqs.  (70-72). 
For  example,  if  the  denominators  of  w^  or  w^  vanish,  Eq.  (72) 
is  appropriately  replaced  by 


JMAX  - 1 

<t)  OMAX  - (j-l) 


and 


= IMAX  - i 
Z IMAX  - (i-1) 


(73) 
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SECTION  V 


CONCLUDING  REMARKS 

In  summary,  the  present  study  represents  a closer  analysis 
of  the  numerical  coordinate- transformation  technique  and  pro- 
vides important  refinements  in  the  procedures  involved  in 
determining  surface-oriented  coordinates  for  arbitrary  geome- 
trical configurations . These  refinements  lead  to  a more 
desirable  distribution  of  mesh  points  for  the  finite-difference 
solution  of  the  corresponding  fluid  dynamics  problem.  The 
resulting  procedures  are  general,  are  applicable  for  internal 
as  well  as  external  flow  problems  and  become  increasingly 
useful  for  high-Re  flow  calculations.  Thus,  the  study  leads 
to  an  efficient  method  for  generating  suitable  surface-oriented 
coordinates  for  arbitrary  geometries. 

Nevertheless,  possibilities  for  further  improvement  con- 
tinue to  present  themselves.  For  example,  even  though  the 
coordinate  solutions  converge  quite  rapidly  for  the  configurations 
tested,  increased  computational  efficiency  is  possible  with  the 
use  of  variable  time-steps  in  the  ADI  method  of  solution  for 
the  coordinate  equations  (4-7) . Finally,  attention  is  drawn  to 
a concern  that  arises  in  the  use  of  a numerical  coordinate 
transformation  in  the  analysis  of  turbulent  flows.  Since  a 
numerical  mesh-generation  procedure  represents  the  body  surface 
by  segments  of  straight  lines,  the  sharp  corners  formed  by 
these  line  segments  may  lead  to  a distortion,  in  the  resulting 
turbulent  flow  solution  as  compared  to  the  true  turbulent  flow. 
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numerically  transform  the  multi-sided  polygon  representation 

17 

of  the  body  surface  into  a sxngle  straight  line  . Since  the 
Schwarz-Christoffel  transformation  i 

mation,  the  resulting  coordinates  may  not,  in  general, 
provide  appropriate  resolution  of  certain  critical  regions 
of  the  flow.  Use  of  a further  numerical  transformation  may 
be  necessary  to  achieve  the  desired  resolution. 

Finally,  it  is  important  to  note  that  the  present  study 
has  led  to  significant  improvements  in  the  convergence  rate  of 
the  numerical  solutions  of  the  coordinate-transformation 
equations.  Most  configurations  tested  (Fig.  1)  required  about 
20  seconds  on  a CDC  6600  computer  for  an  (11  x 41)  grid  to 
satisfy  a convergence  criterion  of  less  than  one  percent  change 
in  the  relative  value  of  the  solution  at  each  grid  point,  per 
unit  time.  The  maximum  time  required  for  any  case  tested  in 
this  study  is  about  80  seconds  for  obtaining  a satisfactory 
coordinate  system.  These  computing  times  represent  an  order- 
of-magnitude  reduction  in  the  computer  time  required  prior  to 
the  present  refinements. 

Optimization  of  the  grid-point  distribution  is  necessary 
in  order  to  reduce  truncation  errors  in  the  numerical  solution, 
especially  in  region,'.;  where  the  flow  variables  undergo  large 
gradients.  For  complex  flows,  these  regions  are  not  known 
a priori.  Adjustment  of  the  grid  with  the  development  of  the 
flow  field  is  then  necessary.  A self-adjusting  grid  generating 
pirocedure  is  therefore  a desirable  feature  for  further  study. 
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APPENDIX  A 


LOCAL  ORTHOGONALITY  AND  LOCAL  CONFORMALITY 


This  appendix  outlines  the  implications  of  settiiij  to 
zero  the  coefficient  b of  the  mixed- derivative  terms  in  the 
coordinate  equations  ( 4-7),  i.e., 

This  condition  can  be  rearranged  as  follows: 


or 


3 <|)  I 3 <t> 

anlj 


3z  3z  I 

an  ^ ac  I 


[a4'/3n]^  [az/ac]^ 

nVaTT^  ==  " ThThT^ 


(A-2) 


Using  the  chain-rule  for  differentiation  enables  rewriting 
Eq.  (A-2)  as 


- - 


(A- 3) 


This  clearly  implies  that  the  curve  c=constant  is  orthogonal  to 
the  curve  n=constant.  Therefore,  Eq.  (A-1)  implies  local 
orthogonality  of  the  (n,C)  coordinates. 

Next,  it  is  shown  that  the  condition  for  the  (n,?)  coordinates 
to  be  conformal  is  a special  case  of  Eq.  (A~l) . Since  a mapping 
defined  by  an  analytic  function  is  conformal,  z+i(})  and  C+in  are 
considered  to  be  complex  variables  related  by  an  analytic  trans- 
formation as 
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z + i(J) 


f(C+in) 


(A- 4) 


Then,  by  definition  of  analyticity,  the  real  and  imaginary  parts 
of  the  mapping  function  satisfy  the  Cauchy- Riemann  conditions, 
so  that 

z,  = (J)  and  z = - 4,  . (A-5) 

C n n C 

Satisfaction  of  these  conditions  implies  the  satisfaction  of 
condition  (A-1) , but  the  converse  is  not  true.  Conditions  (A-5) 
place  greater  restriction  on  the  transformation  than  does 
condition  (A-1) . This  can  also  be  seen  by  observing  that  conditions 
(A-5)  imply,  further,  that  ij>  and  z are  harmonic  functions 
satisfying  the  Laplace  equation.  In  fact,  use  of  conditii  ns 
(A-5)  in  Eqs . (6-7)  leads  to 

a = c = J and  b = 0 . (A-6) 

Furthermore,  the  inverse  transformation  corresponding  to  Eq.  (A-4) 
is  also  analytic  and  conformal.  Therefore,  c and  n satisfy  the 
Cauchy- Riemann  equations. 

and  ^ ~ 

Hence,  n and  ? must  be  harmonic  functions,  so  that  the  forcing 
functions  Q and  R in  Eqs.  (1-2)  must  be  identically  zero.  With 
zero  Q and  R and  with  Eqs.  (A-6) , the  inverted  coordinate  equations 
(4-5)  reduce  to  simple  Laplace  equations.  This  is  clearly  more 
restrictive  than  desired  for  a general  configuration  where 
computational  efficiency  and  accuracy  may  require  a highly  non- 
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uniform  mesh-point  distribution  in  the  physical  plane  and, 
hence,  non-zero  Q and  R in  the  numerical  coordinate-transfor- 
mation. 


APPENDIX  B 


INVERSION  OF  EQUATIONS  (42a,b) 

This  appendix  presents  the  details  of  the  derivation  of 
Eqs.  (46-47)  which  represent  the  inverted  forms  of  the  equations 


^N,S  “ Q 


"n,s  ^ ^ 


where 


N,S 


1 ''a  '^a^ 

1[N  1-J+  H |jT+S 

h^  3N^  as^ 


' ' a 

" ® Is^ 


(B-1) 


(B-2) 


These  are  Eqs.  (42-43)  in  the  text.  The  symbols  used  in  this 
appendix  have  been  defined  in  Section  3.2.1  of  the  text. 
Starting  with  the  functional  relationships 


N = N(n,?) 


and 


S = S(n,c) 


{B-3) 


for  the  numerical  coordinate-transformation,  the  various  derivatives 

2 

appearing  in  the  operator  _ can  be  expressed  in  terms  of  their 


inverted 

forms , as 

follows ; 

^N  = 

’ s^/j  , 

= -S  /j  , 

N n ' 

-N^/j  . 

's  = "/j  ' 

(B-4) 

’^NN 

- - 

S^A2)/J 

''NN 

+ , 

^SS 

= - 

S^B2)/J 

'“SS 

+ , 

(B-5) 

where 


^^2  = Wn.^c  ^ \c®n  - 


B,  = (S  N?  + S_,N^ 

1 nn  ? ^ <,  A 

B.  = (N  N?  + N^  N^ 

2 nn  C C?  n 


2N^^N^N^)/j2 


(B-6) 


and 


J = N S,  - NS 
ri  C C n 


(B-7) 


The  derivation  of  Eqs . (B-4  - B-7)  parallels  that  described  in 

Refs.  4 and  5.  Substitution  of  these  equations  into  Eqs.  (B-1  - 
B-2)  leads  to  the  two  equations  given  below; 

2 2 

n'  (N^A^  - S^A^)  + n"s^  + s'  (N^B^  - S^B2)  - s"Nj.  = J h^Q 


and 

.2 

N (-N^A,  + S^A.) 
n 1 n 2 


2 

n''s  + s'  {-N  B,  + S B-)+  s"n  = j h^R  . 

n n 1 n 2 n 


(B-8) 


These  are  rearranged  to  yield  the  following  form; 


2 2 2 1 2 , 

(n'  + s'  B^  - s")N^  - (n'  A2  + S B2  - n”)Sj,  = j h^Q 

and 


^ « 2 2 

-(n'  At  + s'  B,  - s'')N  +(N  A^  + S B„  - n'')S  = J h^R  (B-9) 

1 1 ' n 2 2 n 


Equations  (B-9)  can  be  viewed  as  two  independent  equations  for  the 

two  quantities  in  parentheses o Therefore,  Eqs.  (B-9)  yield 
2 2 

n'  a,  + s'  B„  - n"  = j h^  (QN  + RN  )/(N  S - NS)  (B-10) 

4^  z n > ^ n M ^ 
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NA,  + SB,  - S =J  h^(QS^  + RSJ/(N.S^  - N„S_)  . (3-11) 

1 1 nCCnnC 


Use  of  Eqs.  (B-6  - B-7)  in  Eqs . (B-10  - B-11)  leads  to  the  inverted 
equations  corresponding  to  Eqs.  (B-1  - B-2)  as  follows; 


2 2 

N*  - 2N^  S„SJ+  s'  - 2N  _N^NJ 

nn  C n nc  n c nn  c CC  n nt  n C 


- n''  + h^(QN^  + RN^)  = 0 


(B-12) 


and 

.2 


.2 


N (S  S?  + - 2S  S S^)+  s'  (S  N?  + S^  - 2S  N ) 

nn  5 55  n nC  n C nn  c CC  n n?  ti  ?' 


- s''  + ° 


(B-13) 


The  inverted  equations  (46a, b)  in  the  text  are  simply  a rearranged 
form  of  Eqs.  (B-12  - B-13)  as 

a N + 2b  N + c N,,  - n"  + h^  (QN  + RN^) 
nn  n?;  cc  n c 

and 

as  + 2b  S , + c S,  - s"  + h^  J^(QS  + RS J 
nn  nc  cc  n c 

where  a,  b,  c are  as  defined  in  Eq.  (47)  in  the  text. 


= 0 (B-14) 


= 0 (B-15) 
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